Dynamics of granular avalanches caused by local perturbations 
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Surface flow of granular material is investigated within a continuum approach in two dimensions. 
The dynamics is described by a non-linear coupling between the two 'states' of the granular material: 
. a mobile layer and a static bed. Following previous studies, we use mass and momentum conservation 

to derive St-Venant like equations for the evolution of the thickness R of the mobile layer and the 
profile Z of the static bed. This approach allows the rheology in the flowing layer to be specified 
independently, and we consider in details the two following models: a constant plug flow and a 
linear velocity profile. We study and compare these models for non- stationary avalanches triggered 
by a localized amount of mobile grains on a static bed of constant slope. We solve analytically the 
^ non-linear dynamical equations by the method of characteristics. This enables us to investigate the 

temporal evolution of the avalanche size, amplitude and shape as a function of model parameters and 

S initial conditions. In particular, we can compute their large time behavior as well as the condition 
for the formation of shocks. 

PACS numbers: 45.70.Ht — Avalanches, 45.70.-n — Granular systems, 05. 45. -a — Nonlinear dynamics and 
£h ' nonlinear dynamical systems 

o 

I. INTRODUCTION 

i ! 

The dynamics of granular avalanches has been keeping busy a large fraction of the granular community for several 
years. As an example, a recent review paper has been written by the French research group named 'GdR Milieux 
' Divises' in order to sum up the 'French results' of experiments and numerical simulations on steady uniform dense 
granular flows . Some of the recurring issues addressed in this long article concern the definition and the description 
of the rheology and the friction law of these flows. As a matter of fact, these quantities are important physical 
ingredients to be plugged into the equations proposed for the modeling of these avalanches. Although some alternative 
models are proposed - see e.g. - an interesting theoretical framework for the description of granular flows follows 

the St Venant-like approach for thin flows in hydrodynamics, in which conservation equations are integrated over the 
depth of the flow Q. This was proposed for example by Savage, Hutter and co-workers for the case of the motion of 
grains over an inclined rough plane p|. 

With such equations adjusted on uniform steady flows, it is for instance possible to reproduce quantitatively steady 
fronts of granular avalanches on a rough inclined plane 6] . Another interesting situation is that of unsteady avalanches, 
whose dynamics is of course more complicated to capture and is very demanding for the models. An example of such 
O ■ a situation is the motion of an initially confined granular mass which is released and runs down the same rough plane. 
In this case, the same calibration gives also quite good qualitative predictions 0- The spreading of a granular mass 
is also of interest in a geophysical context as they are a model of real cliff collapse 0, 0] . Another famous example is 
the case where an initial static layer of grains is available on the plane, and avalanches triggered by a local solicitation 
with a pointy object |ld lllj| . The nice experimental pictures show two types of behavior: a 'triangular avalanche' 
when the thickness of the grain layer is small, and an avalanche with an 'up-hill' propagative front when the layer is 
thicker. 

In the present paper, our aim is also to deal with unsteady and non-uniform granular flows in a geometry close to 
these elementary 'response' situations where avalanches are caused by local perturbations. More precisely, we look 
at the dynamical evolution of a uniform slope on which some amount of rolling grains is initially deposited locally. 
However, unlike the above examples, we shall focus on avalanches running over an erodible bed, i.e. an infinite grain 
layer for which the selected thickness of moving grains is not fixed a priori or bounded by a plane. This situation 
then has one extra dynamical degree of freedom. As will be made more precise in the next section, in order to close 
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FIG. 1: Rolling grains are locally deposited on a static sand bed which has initially a uniform slope p + m, where p is the 
tangent of the angle of repose. 



these equations it is indeed necessary to specify the shape of the velocity profile of the moving grains and to explicit 
the different forces acting on them. Moreover, as our theoretical work aims at obtaining exact results on avalanche 
profiles, the analysis will be restricted to two dimensional systems for simplicity, and will be relevant for experiments 
in confined cells or when the perturbation has translation invariance in the direction perpendicular to the slope. 

General St Vcnant equations for granular avalanches on erodible beds have been recently proposed by Douady et 
al. in [l^. Observations show that these avalanches consist of a thin moving layer of grains over a (quasi)static profile 
|l3j. The two natural variables in this modeling are therefore the profile of the static pile Z, and the thickness of 
rolling grain layer R that flow over this bed, both functions of time and horizontal position. In this framework, we 
shall derive two particular sets of coupled equations for R and Z that are simple enough to be solved exactly for the 
'response geometry' described above. As will be explicited in the next section, they are based on different choices 
for the velocity profile in the flowing grain layer, and exhibit different non-linearities. A particular case of one of the 
two studied sets gives in fact the equations proposed in on more phenomenological grounds, and for which we 
developed an analytical treatment in 0] using the method of characteristics. We extend here this powerful technique 
to all other cases as well. This technique has been also applied to a St Venant like model for debris avalanches [lq . 

After the derivation of these sets of equations in section [D] in which we shall specify the assumptions made, both 
sets of equations - corresponding to different flow profiles - are treated in the following two sections. We show the 
different predictions of the models, compare how the triggered avalanche dies out or grows depending whether the 
initial slope of the static profile is smaller or larger than the repose angle of the pile, and discuss the appearance 
of shocks. We finally conclude and draw perspectives on the description of avalanche fronts. The present paper is 
rather technical, but we hope that the scope of the method used is sufficiently broad to warrant interest in itself; 
furthermore, obtaining explicit exact result for the shape of avalanche shapes and sizes gives important benchmarks 
to which experiments can be compared, and will eventually help selecting the correct set of hydrodynamical equations 
for granular avalanches. 



II. ST VENANT EQUATIONS FOR GRANULAR AVALANCHES 

The aim of this section is to (re-)establish coupled differential equations for the variables R and Z introduced above 
and depicted in figure [3 These equations encode the conservation of the number of grains as well as their horizontal 
momentum |l2j . They describe how static grains may be dislodged and contribute to the flow (erosion) , and vice- versa 
how moving grains may come to rest (deposition) |14j . 

Besides R and Z, another important field is of course the velocity in the moving layer. Let us note u(x, y) the 
horizontal component of the velocity profile in this layer - although not explicited, time dependence of u is understood. 
The averaged velocity of the flow can be defined as U = J^ +R dyu(x, y). With this quantity, the number of particles 
passing through a vertical line in x during the time interval dt is ^RUdt, where p is the mass density of the granular 
material and m the mass of a single grain. The difference of this number in x and x + dx makes the volume (R + Z)dx 
change, so that the conservation of matter reads 

d t (R + Z) = d x (RU). (1) 

In a similar way, the ^-component of the momentum passing through the vertical line in x during dt is pRWdt, where 
W is the average of the square of the velocity: W — J^ +R dy u 2 (x,y). The change of horizontal momentum pRUdx 
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FIG. 2: Rolling phase R(t,x) (grey) and static phase Z(t,x) (crosses). The grains flow from the right to the left. The dashed 
box is the control volume for which mass and momentum conservation laws are written. 



in the control volume (see figure^ can then be written as 

d t (RU) dx = d x (RW) dx + - dF x , (2) 

P 

where dF x is the ir— component of the forces acting of this volume. 

These forces are of two kinds: those coming from the lateral sides through the stress a xx , and those due to the 
contact between the rolling and the static phase. Assuming a horizontal normal stress proportional to the 'hydrostatic' 
pressure o xx (x, y) = b pg(R + Z — y), the corresponding force is J^ +R dya xx = \bpgR 2 . Numerical simulations suggest 
values for b that increase from roughly one half in static piling to unity for uniform steady flows 0, ITtI ITsj . In [j| the 
value of b is obtain for unsteady and inhomogeneous situations from a Mohr-Coulomb yield criterion. It turns out 
that b can increase (decrease) compared to unity if the grains are compressed (decompressed) by the flow. Here we 
consider b as a phenomenological constant. Taking the difference of this force on the two sides of the control volume, 
the contribution of the horizontal normal stress to dF x finally reads bpgRd x Rdx. To compute the other part of dF Xl 
we assume that the grains in the control volume behave, with respect to the static phase, like a frictional solid (with 
normal reaction N and friction force T). We call p the friction coefficient. Defining 9 by tan# = d x Z, the balance of 
the weight of this volume gives pgRdx — N cos 9 + T sin 9. With T = pN, the horizontal contribution of N and T to 
dF x (to the lowest order in 9 and p) is pgR{d x Z — p) dx. In summary, we get 

dF x = bpgRd x R dx + pgR(d x Z — p) dx. (3) 

Several remarks are in order at this point. First of all, we have assumed that the granular density p is the same 
in the moving and the static parts, which is certainly a valid hypothesis for dense flows as considered here. More 
importantly, the friction coefficient p will be taken as a constant in the following. It is however well established that 
there are few degrees of hysteresis between the starting and the stopping angles, which can be simply understood even 
at the scale of a single grain rolling down a "pile" consisting of a layer of regularly spaced fixed grains 0] . As will be 
emphasized in the conclusion, this slight difference is of great importance when it comes to the description of regions 
where the moving layer is about to jam or the static bed about to move, i.e. at the feet of avalanche fronts (R — > 0). 
Besides, it has been shown that in quasi 2D experiments conducted in a thin channel between two plates, the friction 
of the rolling on those plates plays a major role in the dynamics of the avalanche |20|,|2lJ. Such a boundary effect is 
not encoded here. 

Finally, these equations will be closed by making an explicit choice for the velocity profile u. The integrals for the 
computation of U and W will then be expressed as a function of R. This means that we implicitly assume instantaneous 
adaptation of the velocity to the flow. In principle, like R and Z, the averaged velocity U is an independent dynamical 
variable whose dynamics should be described by an additional equation of 'internal force balance' dtU — .... The 
right hand side would specify the part of the forces acting on the moving layer that contribute to the acceleration and 
friction of the rolling grains rather than the erosion/deposition processes, i.e., the exchange between R and Z . A model 
similar to our Eqs. CEJ - © has been very recently studied by Khakhar et al. but for different geometries and 
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initial conditions (heap and rotating cylinder flows). Their momentum balance equation has however an additional 
term which describes grain collisions in the mobile phase. Thus they combine effects of internal dissipation and 
external forces as gravity and solid like friction between the two phases. 

Two choices for the velocity profile will be investigated below: the plug flow for which u(x, y) = U = Cst (named 
model V for 'plug'), and the situation with a linear velocity profile u(x,y) — jy (named model C for 'linear'). The 
natural selection of this profile is still a puzzling issue and is the result of several mechanisms at the grain level (e.g. 
trapping) or at larger scales (non-local effects, clustering) At the phenomenological level, the first case is 

reasonable for the description of thin or dilute flows (when R of the order of a grain diameter) [26| , and the second one 
is strongly supported by experiments and simulations performed on steady and deep 2D systems pj. In particular, it 
should be emphasized that the velocity gradient is found to be constant, and it is the thickness of the rolling layer R 
which adapts its value to the external imposed shear stress, see e.g. [l|,|27j. In contrast, note that for steady granular 
flows on a fixed rough inclined plane a Bagnold-like velocity profile U cx i? 3 / 2 is observed. For the study of unsteady 
situations, since the static profile Z is now specified, the two St Venant conservation equations can be used for the 
determination of the space and time evolution of R and U Q. 



For a plug flow with a constant average horizontal velocity U we trivially get W = U 2 , and one obtains from the 
analysis in section [H] the model 



Here the dimensionless t, x, R and Z are measured in units of U 2 /g, and t is rescaled by U/g. Furthermore, as we 
are interested in surface profiles close to the avalanche slope /i, Z is measured relative to the critical slope, i.e., it 
is replaced by Z + fix. The quantity b is then the only free parameter of these equations. Recall that an isotropic 
stress distribution corresponds to b = 1. Note also that setting b = yields the so-called BCRE model introduced 
in (although without the diffusive terms considered there). Any small but finite b thus leads to novel non- 
linearities. Before we study the propagation of a localized perturbation, as a first simple benchmark of the model, we 
briefly note the predictions of the model for the initial situation of a constant slope Zq(x) = mx and a homogeneous 
amount of rolling grains, Rq(x) = g. It is easy to see that in this case the thickness of the mobile layer growths (or 
decays) exponentially in time, depending on the sign of m, with R(t,x) — ge mt . The static profile is then given by 
Z(t, x) — mx + g(l — e mt ). Note that this solution does not depend on b since R(t, x) is independent of x. 



This non-linear model has been studied previously for the case of a static profile consisting of two regions of constant 
but different slopes and an initially homogeneous (constant) amount of rolling grains [15|. However, the method of 
characteristic curves can be also employed to study local perturbations in form of an initially localized amount Rq(x) 
of rolling grains. In fact, for arbitrary initial profiles Rq(x) and Zq(x) the solution of equations I0J with 6 = can 
be obtained analytically in implicit form. Using the method of characteristic curves, see Appendix [XJ we introduce 
the new coordinates a(t,x) and j3(t,x). The mapping between the two coordinate systems depends on the solution 
for R(t,x), Z(t,x) and thus on the initial profiles due to the non-linearities. The so-called characteristic equations, 
which determine simultaneously the coordinate mapping and the profiles, for this model read 



III. PLUG FLOW: CONSTANT VELOCITY PROFILE (MODEL T) 



d t Z = -Rd x Z -bRd x R, 
d t R = d x R + Rd x Z + bRd x R. 



(4a) 
(4b) 



A. Infinite stress anisotropy (b — 0) 



1. Analytical solution 



d a x — (,+d a t = 0, 

dpx-Q-dpt = 0, 

Rd a Z + (C+ - R) d a R = 0, 

RdpZ + (C_ - R) dpR = 0, 



(5a) 
(5b) 
(5c) 
(5d) 



with the characteristic directions given by 



C+ = -i, C-=R- 
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FIG. 3: Bold curve: plot of the Lambert's function W(x). Thin curves: small and large x expansions W(x) ~ x — x' 2 and 
W(x) ~ In a; — lnlnx, respectively. 



The solution of these equations can be obtained explicitly, 



dp' 



t(a,f3) = - 

J -a 1 + R{<X,P) 

x(a,/3) = -p-t(a,0). 
As function of these new coordinates the solution can be written in the closed form 

Z(a,/3) = Z (a) 



(7a) 
(7b) 

(8a) 
(8b) 



where W is Lambert's function (2^, see Fig- EI Before we proceed with special choices for localized initial profiles 
Rq(x), we would like to point out that the model of equations J3J with 6 = can be still solved exactly if the coefficient 
of d x Z on the right hand side of both equations is a general function of R(t,x), T{R). The above case corresponds 
to T{R) = R. 

In order to progress with analytical techniques we make a special choice for a localized R${x) which allows for a 
closed expression for the integral of the coordinate map in equation Ij7a(l . Below, we will demonstrate by an explicit 
numerical integration of the coordinate map for a generic Gaussian perturbation Rq{x) that the following results are 
robust with respect to the precise form of the perturbation. Thus we proceed with the choice 



Ro(x) = W 



r e 



ro-\x\/S 



(9) 



for the initial profile. This profile decays exponentially at large \x\ and has a amplitude of ro- For the relevant limit 
of ro <C 1 the width at half amplitude is 6 In 2. This form is adapted to the integral in equation (|7a|l since the latter 
equation can be written as 



t(a, 0) = —a — (3 + 
and Rq (x) of equation has the property that 



df),R(a,P')dP' 



a cV(lni? (-/?')) - R (-(3>) - Z Q (-f3>) 



(10) 



d x (lnRo(-x)) - R' (-x) = --sgn(x). 

o 



(11) 



For the static phase we will consider always a profile with a constant slope, 



Zq(x) = rax, 



(12) 
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where due to our definition of Z the parameter m measures the excess slope relative to the critical angle. For these 
initial data the solution of equations Q in the curved coordinate frame reads 



R(a,/3) = W 



r e 



ro-\0\/6-m(a+P) 



ma. 



(13a) 
(13b) 



In order to integrate the equations for the coordinate mapping, see equation Q7a[l. we have to divide the space-time 
into different sectors due to the sign function of equation (|llfl . Depending on the sign of a and we define the 
following three sectors, (I): a < 0, (3 > 0, (II): a < 0, (3 < 0, and (III): a > 0, (3 < 0. The case a, (3 > is mapped 
to negative times, and is therefore not of interest. The boundary between the sectors (II) and (III) in the t — x plane 
can be obtained by integrating equation H7a|) with a = and /3 < 0. The boundaries are given by 



xyuit) 
Xu/uxit) 



'o 



,{m-l/S)t 



?71—l/(5 V 

where the subscript indicates the adjacent sectors, see figure 0J 



(14a) 
(14b) 




FIG. 4: Model V: Boundaries between different sectors in the t — x plane. 



In sectors (I) and (III) an explicit expression for the profiles can be obtained. From the integrated version of 
equation l|7a|) one easily gets the result 



R(t,x) = R {a{t,x))e (m±1/S)t , 



(15) 



where ± refers to sector (I), (III), respectively. Together with equation (|13b|l this result shows that the system has a 
simple time evolution along the characteristic curves of constant a(t, x). Using equations Ijl3bf> . (|15f) and (3 = — t — x, 
we obtain the explicit expression for the coordinate mapping in sector (I), 



— Z(t, x) = a(t, x) = x + 
m 



D {rn+l/5)t 



m+ (l/<5)e( m + 1 /<5)* 



w 



ro 



m +1/5 



e ro+x/5 m . 



(l/5)e 



(m+l/S)t 



(16) 



The result for sector (III) is obtained from the latter expression by the replacement S — * —S. equations (|15l) . Ijltjfl are 
our final result for the profiles in sectors (I) and (III). In sector (II) the characteristic curves x a (t) which maps to a 
constant a can be obtained again explicitly but we were unable to invert them to obtain a(t,x). The characteristic 
curves read 



x a {t) 



m- 1/8 



/i(a)e (m " 1/,5)t + \n(h(a)/r ) + ma - r 



with the function 



h{a) = W^T* [r e r °- ma ] 



r e 



r +a/S 



(17) 



(18) 



As in sectors (I) and (III) the curves behave exponentially in time with, however, more complicated amplitudes. The 
characteristic curves are shown for ro = 0.1, S = 5.0 and different slopes m in figure |3 

With the characteristic curves at hand, the solutions for i?(t, x) and Z(t, x), at a fixed time, are given in parametric 
form in the x — R or x — Z plane by the curves [x a (t), R(a, (3 = —x — t)], [x a (t), ma], respectively, with a as running 
parameter. 





FIG. 5: Model V: Characteristic curves with constant a for a perturbation Ro(x) of Eq. @ with ro = 0.1, 8 = 5.0 and slopes 
(a) m = —0.05, (b) m = 0.05 < 1/8, (c) m = 0.15 < 1/8 and (d) m = 0.25 > 1/8. The shaded area indicates the region with 
shocks. The dashed curves mark the boundaries between the different sectors, cf figure 0] 



2. Appearance of shocks 



Before we discuss the resulting profiles, we will study the possibility for the occurrence of shocks, i.e., discontinuities 
of the profiles which develop at a finite time. Such kind of singularities are possible for non-linear dynamics since 
adjacent characteristic curves can bend differently, leading even for small differences in curvature at small times to 
intersecting curves at larger times. Beyond the time at which the curves cross for the first time, there is a region 
where a unique solution no longer exists. Two such situations are visualized for our model in figures [5f a) and (c). If 
characteristic curves cross each other they form an envelope. The shock position is then determined by the cusp of 
the envelope. While we leave the precise definition and the calculation of the envelope to Appendix^! we discus here 
the criterion for the existence of shocks in our model. We are interested in not too large amplitudes ro. For ro < y/2, 
there is a simple condition for the formation of shocks. Then it can be shown, see Appendix 151 that shocks occur 
only if the slope of the profile of the resting grains is larger than the critical value 

m c = ^fl. (19) 



Thus for broad perturbations in the moving phase shocks occur already for small slopes to. The exact time t s 
and position of the shock is in general difficult to calculate. For sufficiently large to (to > l/(2ro<5) if r < 1 and 
to > 1/(5(1 + ro)) if ro > 1) an exact expression can be obtained and is given in Appendix iBl On the other hand, for 
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slopes close the critical one the shock time diverges as 

t s ~ (m-nic)" 1 . (20) 
For our above choice of parameters ro and 5 for figure[S]the critical slope is m c = 0.0828. 

3. Different avalanche dynamics 

Figure [S] shows four possible situations: (a) A slope to which is larger than critical m c but smaller than 1/5 so that 
the characteristic curves saturate at large times, (b) a slope smaller than the critical m c , (c) a slope which is larger 
than both the critical m c and 1/5 so that the curves grow exponentially in time, and (d) the case of a negative m. 
The solutions for the corresponding profiles of R(t,x) and Z(t,x) are shown in parts (a) of figures El - |H1 Plotted 
are the two curves Z(t,x) — mx (always corresponding to the lower curve) and Z(t,x) — mx + R(t,x) so that the 
gap between the two curves represents the layer of moving grains. The maximum of the perturbation propagates 
downhill with a constant velocity which is 1 in our rescaled units. For a negative slope to (corresponding to an actual 
slope smaller than the angle of repose) all grains of the perturbation come finally to rest, generating a bump on the 
initial static profile which corresponds to the baseline in the plot, see figure Efa). The amplitude of the perturbation 
decays exponentially at large times, R(t,x = —t) ~ roexp(ro + rat). If we define the downhill end of the bump by 
the condition that the maximum of the perturbation in R{t, x) has decayed to some fraction e < 1 of the initial 
amplitude, then the width of the bump scales like m(e)/m for small r since the peak in R(t, x) moves with a velocity 
of one. Thus the final width of the deposited amount of grains is independent of the amplitude and width of the 
perturbation but only determined by the initial slope of the static phase. For positive slopes to the amplitude i? ma x 
of the perturbation grows at large times linear, 

i? max ~ m(l + mS) t, (21) 

with a growth rate which is independent of the initial amplitude ro. However, there is a broad transient behavior 
with logarithmic corrections at intermediate time scales. The exponential time behavior found at the beginning of 
this section for an initially homogeneous mobile phase is recovered for m < but is in contrast to the linear growth 
for positive m. 

For both positive and negative initial slopes, the profile Z(t, x) of the static phase does no longer evolve after the 
avalanche has passed. The resulting profile Z^ix) becomes thus time independent for times larger than a (position 
dependent) time scale. This asymptotic profile is only well defined if no shock occurs. In the latter case it is implicitly 
given by 



1 



m — 1/5 



In ( ) + Z 00 (x) - r 



(22) 



where the function h is given by equation l|18|l . The latter expression is valid in sector (II) which is the relevant region 
for large times, see figure From this result one can obtain the slope c^Zqo at a given value of Z, 



dxZooix) = (m - 1/5) 



h'(Z OQ (x)/m) ^ 
mh{Z 00 {x)/m) 



(23) 



The behavior of 2 M depends on the sign of to. For negative to the deposited grains form a bump that was described 
above and whose exact shape is given by equation Q22[l. Across this bump, we find that the change of the slope as 
compared to the initial slope to is always rather small. Far away from the bump (at large positive and negative x) 
one has h'/h ~ —1/5, and thus d^Z^ = m remains, of course, unchanged from the initial profile. For positive to 
the profile Z(t,x) shows again a constant slope after the avalanche has passed, i.e., = m^x, cf figure |7| The 
asymptotic slope can be obtained from the behavior of the function h at large negative arguments. We find 



1 



2m 



i5 + 1 + a/2 



(24) 



Interestingly, the relative change of the slope is independent of the amplitude ro of the perturbation and depends only 
on the product of the initial slope to and the width 5 of the perturbation. Surprisingly, the expression in the square 
brackets in equation is larger than one for positive to, and even diverges as the slope to approaches the critical 
value to c beyond which shocks occur, TOoo ~ (m c — to) -1 . It is important to note that the layer of moving grains 
decays to zero at large times at the uphill end of the avalanche, although the slope of the static layer is steeper than 
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FIG. 6: Model V with b = 0: Fixed time profiles Z(t, x) — mx and R(t, x) + Z(t, x) — rax so that the gap between two 
corresponding curves represents the thickness of the mobile phase. Plot (a) is for Ro(x) of Eq. © and (b) for a Gaussian 
Ro(x), both for m = —0.05. The profiles are for times t = 0, 5, 10, 20 and 50. Here and in the following plots of profiles the 
plotting intensity decreases with increasing time. 

before the avalanche started. This behavior can be easily understood from the property of equation 0} in which the 
exchange between the static and rolling phase is proportional to R. Of course, physically, there will be a maximum 
angle for the stability (even for R — 0) beyond which the above predictions are irrelevant, and an extended model 
has to be considered. 

An important quantity is the total size of the avalanche. Within our model, we define the size I(t) of an avalanche 
as the spatially integrated amount of mobile grains, i.e., 

/oc 
R(t,x)dx. (25) 
-oo 

The integration can be performed for each sector separately by a change of variables to the characteristic coordinate 
a. The contribution from sector (III) can be neglected at larger times since the avalanche starts at x = to propagate 
to the left (downhill), cf figure 01 For negative to, the perturbation decays exponentially, and thus we obtain for the 
size 

J(t) = 2r Se ro+mt . (26) 
For positive to we observe that the size of the avalanche shows a quadratic increase in time, 

m =mH l -±^-t\ (27) 

1 — TOO 

at asymptotically large times. Interestingly, the growth of the avalanche depends only on its initial width 8 but not 
on the amplitude r . By comparison with the scaling of the amplitude of the avalanche, cf equation (|21|) . we observe 
that the width of the avalanche must grow also linear in time ~ mS/(l — mS)t. 

So far we have studied mainly an initial perturbation which is particularly suited for obtaining analytical results. In 
order to check the robustness of our results with respect to the precise form of Rq(x) we have chosen also a Gaussian 
Rq(x) = ro exp(— x 2 /S 2 ) together with the same Zq(x) = mx as before. Contrarily to the previous case, the initial 
perturbation has no cusp at x = 0. By a numerical computation of the integral of equation (|7all we obtained the 
profiles shown in parts (b) of figures |S] - |S1 using ro = 0.1, 6 — 5 as before. As can be observed from the plots the 
characteristic features can be regarded as robust. However, the moving layer, i.e., the gap between the upper and 
lower graph decays more rapidly due to the faster decay of the Gaussian profile. Of course, the critical slope m c for 
shocks is no longer given by equation (|19fl . We observe that for the Gaussian Ro(x) shocks occur only beyond a m c 
which is increased compared to the exponentially decaying profile of equation (|5J with the same width at half height, 
cf figure |S| 
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FIG. 7: Analog of Fig. ©but for m = 0.05 and times t = 0, 5, 20, 40 and 60. 
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FIG. 8: Analog of Fig. ©but for m = 0.15 and times t = 0, 10, 20, 40 and 60 for (b) only. In plot (a) a shock occurs at the 
uphill end of the avalanche (shock time t s — 41.71). Note that for the Gaussian Ro(x), plot (b), shocks are generated only for 
larger values of m, see text. 



B. Strong stress anisotropy: small b 

1. Analytical result for general slopes m 

So far we have assumed such a strong stress anisotropy that there is no horizontal stress a xx , i.e., b = in the model 
of equation Q . In this section we will study the influence of a small a xx on the avalanche dynamics we found in the 
previous section. Although steady state simulations suggest a value of b close to one ^3,0], it is interesting to study 
the regime of small b in order to compare to the BCRE model. The method of characteristics can be applied of course 
to arbitrary values of 6, yielding a system of equations which we could not solve explicitly so that we had to resort 
to a perturbative treatment. Thus, we consider the terms proportional to b in equation Q as small perturbations of 
the b = solution. This can be done by the following ansatz 



Z = Zx + bZ 2 
R = R\ + bR.2, 



(28a) 
(28b) 



where Z\ and R\ denote the solution for b = of the previous section. Although one expects realistic values of b of 
order unity, the perturbative calculation should allow for a qualitative assessment of the effect of a finite horizontal 
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FIG. 9: Model T with finite b = 0.5: Profiles for a Gaussian R (x) and (a) m = -0.05 (for t = 0, 5, 10, 20 an d 45), (b) 
m = 0.05 (for £ = 0, 10, 20, 30 and 45). The dashed curves represent the corresponding profiles for 6 = of section Hll Al and 
are shown for comparison. 



stress. The dynamics of the corrections are then described by the following linear coupled equations 

t Z 2 = -R 1 (d x Z 2 + d x R 1 )-R 2 d x Z 1 

d t R 2 = d x R 2 + R 2 d x Z 1 + R l {d x Z 2 + d x R 1 ) . 



(29a) 
(29b) 



Since we consider corrections of linear order in 6, the terms containing derivatives of the profiles R 2 and Z 2 have 
exactly the same form as those in equation J1J with 6 = 0. Thus, the characteristic directions C+ = — 1, C- = Ri an d 
characteristic curves remain unchanged for small 6. In terms of the characteristic coordinates a and /?, the corrections 
obey the equations 



d a Z 2 



1 + Ri 



d a R 2 



d x Z 1 + d x R 1 )d a t = 



d Z 2 + {R 2 d x Z 1 + R l d x R l )d t = 0. 



(30a) 
(30b) 



Here all functions have to be considered as depending on a and /3, in particular t(a,/3) is given by equation J7a|. In 
order to express the derivatives with respect to x as functions of a, (3, we use d x = d x ad a + d x /3dp together the 
relations 



d x f) = -l, d x a = 



1 



1 



1 + Ri d a x 



(31) 



with x(a,/3) given by equation l|7bj) . Using these relations and the solutions for b = of equation JBJ with the initial 
condition Zq(ol) — ma, the equation (|30f) can be rewritten after integration over a as 



Y I p 
Z 2 + - 1 R 2 + In 



Ri 



l + Ri 

i + M-P) 



d Z 2 + Y^dpRi 



-0 



da'd Ri(a',f3)d a/ x(a',(3) = 



(1 + Ri) 2 d a x \ l + Ri 



Ri 



i?2 



0, 



(32a) 
(32b) 



where the functions depend on a and (3 unless arguments are written explicitly. The explicit expression for d a x can 
be obtained from the solution for 6 = 0, leading to 



d a x — 



1 



l + R (a) 



p d a Ri(a,P') 
. a P [l + Ri(a,(3'W 



(33) 



The function R 2 can be eliminated from equation (|32b|) by using equation l|32a(l . and the resulting linear ordinary 
differential equation for Z 2 can be integrated easily. The result is 



Z 2 (a,f3) = exp 



d(3'gi(a,[3') 



d/3'g 2 (a, (3') exp 



0' 



d0" gi (a,0") 



(34) 
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with the functions 

mi?i 

91<I> " 5) = -(l + RrfdaX (35a) 



This is our final result for Z2, the profile R2 can be computed now from equation Q32a[l . Using the explicit result 
for Ri of equation (|8bll the multiple integrals can be performed easily numerically. The resulting profiles are shown 
in figures [5] for the parameters ro, S and m of section All Al with b — 0.5. For comparison the solution for 6 = are 
also shown as dashed curves. As can be seen from figure EJa), for negative m the avalanche and the static profile are 
not much effected by the presence of horizontal stress ~ b. In contrast, for positive slopes m there is more visible 
effect of the horizontal stress. This effect grows with increasing time scale, cf figure OJb). As one could naively 
expect, horizontal stress has the tendency to shift the peak in dynamics phase downhill (to the left). This shift will 
be analyzed quantitatively below for the case m — and by comparing the numerical results of figures l§l 1101 the shift 
appears to be independent of m. For the static profile Z(t, x) the horizontal stress leaves the final slope after the 
avalanche almost unchanged but it produces a steeper slope in Z(t, x) where the moving layer R(t, x) has maximal 
thickness. There is also a tendency for the static profile to form a local dell at the peak of the avalanche, see figure 
|§Ib). This observation becomes especially pronounced for m — 0, see figure ITUl 



2. Explicit results at the angle of repose (m = 0) 

It is obvious form the structure of equations (|34l) - (|35b() that major simplifications occur if the initial static profile is 
exactly at the angle of repose, i.e., m = 0. Thus, we can study easily the effect of a horizontal stress (finite 6) in this 
situation. Let us first summarize the results in the absence of horizontal stress (b = 0) when m = 0. The static layer 
stays for all times at the angle of repose, Z\{t, x) = 0, and the initial perturbation Ro(x) in the moving layer simply 
propagates downhill, Ri(t,x) = Ro(t + x). Thus the amount of grains in both the static and the moving layer are 
conserved. The situation is different for finite b. The layers of static and dynamic grains, respectively, have coupled 
dynamics with the corrections from finite b given by 

Z 2 (a,f3) = R (a) - fl (-/3) + In ( l±^B ) (36a) 

V 1 + R {a) J 

P) = 1 K(-i9)t(a,/3) - Z 2 (a,{3)} , (36b) 

1 + K {-p) 
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where the mapping between t, x and a, (3 has for to = the simple form 

a;(a,/3) = -t(a,/3)-/3. (37b) 

These equations provide a closed parametric form of the dynamics for a general initial perturbation Ro(x). The 
resulting time evolution of a Gaussian perturbation Rq(x) = vq exp(— x 2 /5 2 ) is shown in figure ITU1 Compared to 
the absence of horizontal stress, b = 0, there are a number of interesting novel features. The avalanche amplitude 
increases and the maximum is shifted downhill. The layer of static grains displays a bump at the initial position of 
the perturbation and a dell which propagates downhill close to the maximum of the avalanche peak. At large times, 
equation (|37a|l leads to a ~ x which together with j3 = —x — t and equation (|36[) yields an explicit expression for the 
profiles. For a small avalanche amplitude ro, the position of the peak follows x = — t — 6/2 for a Gaussian Rq(x). The 
maximum of R(t,x) grows linearly with time, R max ~ b(rQ/6)t, while in the absence of horizontal stress (b — 0) it 
remains constant. Notice that this linear growth was observed for b = only at angles larger than the repose angle 
(to > 0), cf equation i|21|) . The form of the static profile Z(t,x) can be directly obtained from equation l|3tja|) . There 
are two identical contributions which are shifted relative to each other by t. The first contribution is approximately 
given by b[Ro(x) — ln(l + Rq(x))]. This term represents the bump at the start position of the avalanche. The second 
contribution has x replaced by x + t, corresponding to the dell traveling with the avalanche downhill, cf figure 1101 
Thus, both structures in the static layer are determined by the initial profile of the perturbation in the moving layer. 



IV. LINEAR VELOCITY PROFILE (MODEL C) 



It has been argued that a constant velocity profile, as assumed in the previous Section, is only applicable to thin 
surface flow |2(j. For a thicker layer of rolling grains, the velocity should depend on the amount of mobile grains. 
Experiments and simulations for steady deep systems suggest a linear profile for the average horizontal velocity 
u(x,y) = 72/ of the flow. With this profile we have U — ^7-R, W = \l 2 R 2 = \V 2 and the conservation conditions of 
equations (JTJ and J5J) yield the model 

d t Z = -d x Z-bd x R, (38a) 
d t R = Rd x R + d x Z + bd x R. (38b) 

Here all lengths are divided by g/7 2 , and time is divided by 7. Again, Z is replaced by Z + \ix. The model contains 
after this rescaling only one free parameter, b. It is rather important to note that the latter model is valid only as long 
as R remains positive since we obtained it from equation © after division by R. Thus the actual solution of equation 
(|38JI is given by the maximum of R — and the formal solution for R of equations l|38a|l and l|38b(l . In this Section we 
will study the consequences of a R dependent linear velocity profile both in the absence (6 = 0) and presence (b 7^ 0) 
of horizontal stress. We note that for an initially uniform amount of rolling grains Ro(x) = g and a static sand bed 
with constant slope Zq(x) = rax the solution to Eqs. I|38(l is rather simple. As opposed to the exponential growth for 
model V 1 the thickness of the mobile layer increases here only linear in time, R(t, x) — g + rat and Z(t, x) = m(x — t) 
decreases accordingly. As for model V this solution is independent of b. 



A. Infinite stress anisotropy (b = 0) 

1. Analytical solution 



In the limit of b = the equations l|38a() and l|38b|) are decoupled. Such set of equations has been studied by de 
Gennes et al. to describe a thick flow of granular matter in a bounded geometry |26| . Here we consider this simple 
model in an unrestricted geometry but we allow for general initial profiles Ro(x) and Zq{x). Following again the 
approach outlined in Appendix 1X1 we obtain the characteristic equations 

(39a) 
(39b) 

-d a Z + (C+ - 1) d a R = 0, (39c) 
-d Z + (C_ - l)d R = 0, (39d) 



d a x - 


- (+d a t 


= 0, 


dpx - 


- t-dpt 


= 0, 


(c+- 


l)d a R 


= 0, 


(C-- 


l)d R 


= 0, 
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with the characteristic directions given in the case of b = by 

C+ = l, C- = -R- (40) 

Since one of the characteristic directions is constant, these equations can be integrated in a way which is similar to 
the procedure we used for model with a constant velocity profile in section IIII Al From this calculation one finds 
easily that the general solution of equations (|38(l with 6 = reads 

J_ a 1 + R[a,P') 

x(a,0) = -0 + t(a,0), (41b) 

R(a,(3) = -1 + y/[Ro(a) + l] 2 + 2[Z (/3) - Z (-a)}, (41c) 

Z(a,f3) = Z (-/3). (41d) 

Studying the configurations we studied in section 1111 Al for model V with a constant velocity profile, we choose 
Zq{x) = tyix so that the integral in equation (|41a|) can be computed explicitly. One obtains, using (|41c(l : 

t(a,f3) = — \y/[Ro{a) + l] 2 + 2m(a + (3)- R (a) - 1 
to L J 

= -[R(a,p)-R (a)]. (42) 
to 

Since in the limit b = equation (|38a|l acquires a simple linear form, we have obviously the result 

Z(t,x) = Z (x-t) =m(x-t). (43) 

At sufficiently large times one has a ~ x + mt 2 /2, and thus equation (|42|l shows that the amount of rolling grains is 
given by 

R(t,x) = R (x) + mt with x = x + — t 2 . (44) 
2. Physical discussion 

From the above result, the shape of the perturbation at a given large time would be the same for all initial slopes 
to. From this result, the maximum of the mobile layer R(t,x) travels with a velocity w max = ?"o + mt/2 which, for 
m > 0, increases linear in time, i.e., the perturbation feels a constant acceleration. This has to be compared to the 
constant velocity we found for the model with a constant velocity profile, see section Till Al However, for negative m, 
as explained above, the actual solution is obtained by setting R{t 1 x) to zero in regions where it would be negative 
otherwise. At the time at which R(t,x) becomes zero, the profile Z(t,x) is frozen at its present height. Due to this 
construction the final solution for Z(t, x) will deviate from the simple form of equation l|43|> . For positive to the profile 
R{t, x) as obtained from equation (|42|) is always non-negative. But there is the possibility that shocks actually prevent 
the system to reach the asymptotic time result of equation 1)44(1. In fact, independent of the initial parameters of 
Rq and to, a shock always occurs after a finite time t s , unless the solution R(t,x) of equation (|44|l becomes formally 
negative (for negative to) before the shock can appear. The shock time is given by the general expression 

ts = — Wt (45) 

max rtg [x ) 

Remember that for the plug flow V model (see section ITTT|) . a shock appears only for slopes to which are larger than 
a positive critical slope. For the model discussed here, shocks can even occur at negative slopes to. 

In the following, we will consider again a Gaussian perturbation in the layer of rolling grains, Rq(x) = 
r exp(—x 2 /S 2 ). Then the shock time is t s = ^fej 2 8/ro . Figure ITTI shows the time evolution of this perturba- 
tion for both positive and negative to. Plotted are the profiles Z(t, x) — mx and Z(t, x) — mx + i?(i, x) so that again 
the gap between the profiles corresponds to the layer of rolling grains. For negative to all moving grains have come to 
rest at the time t — —ro/m which is smaller than the shock time scale t s for the parameters used here. For positive 
to there is an uniform increase in the thickness of the layer of rolling grains, see figure UTT b). This increase is linear 
in time, ~ mt, and is unrelated to the amplitude of the local perturbation Rq(x). This apparently unphysical result 
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FIG. 11: Model £ with 6 = 0: Profiles for a Gaussian profile Ro{x) with ro = 0.7, S — 5.0 and slope (a) m = —0.1 (b) m = 0.1, 
both for times t = 0, 3 and 7. The shock time is t a = v/ e/25/ro = 8.32. For negative m the perturbation has stopped at the 
time t = —ro/m — 7.0, i.e., the profile shown for the latter time is the final static profile. Note that for positive m the thickness 
of the moving layer shows an overall linear increase proportional to mt since for a Gaussian profile, Ro(x) is in fact small but 
finite everywhere. For clarity, we indicated explicitly the thickness R of the mobile layer at t = 3. 



can be understood from the structure of equation (|38bl) . Even for a strictly localized initial Ro(x) which is zero 
outside a finite interval, there would be an increase ~ (d x Z)t (for a constant slope) at all positions x, not only there 
where Rq(x) is non-zero. But since we divided the original equations by R to obtain equation i|38[l . R = is a trivial 
solution. The latter solution should be matched with the finite R solution at the front of the avalanche. However, 
by definition, at the front the rolling layer becomes very thin, and a strictly linear velocity profile is certainly an 
oversimplification. Thus, with the model of this Section, the matching of the two solutions at the front is not justified. 
Instead, the dynamical equations should be refined as to describe the phy sical processes close to an avalanche front 
and the thin-to-thick flow crossover (for example along the lines of pfl l29j). This we leave to a future work. 



B. Strong stress anisotropy: small b 

1. Analytical results 

Now we study the influence of finite horizontal stress with a finite but small b. An important consequence of a 
finite b is that now the equations (|38|l become coupled by the stress term. In order to obtain the dynamic response 
to a local perturbation we perturb about the 6 = solution of the previous section. Following the analysis of section 
IIII Bl we make the ansatz 

Z = Z x + bZ 2 (46a) 
R = Ri+bR 2 , (46b) 

where Z\ and R\ denote the solution for b = of the previous section, cf equations (|41c|) and 141d|) . By expansion of 
the equations l|38() in b we obtain the dynamics of the contributions from finite b, 

d t Z 2 = -d x Z 2 -d x R! (47a) 
d t R 2 = {l + R 2 )d x R 1 +R l d x R 2 +d x Z 2 . (47b) 

To this coupled system of equations we can again apply the method of characteristic curves. The resulting character- 
istic equations for the corrections to the profiles are 



d a Z 2 + d x R 1 d a t = 0, 
dpZ 2 + (1 + -Ri) dpR 2 -[(1 + R 1 )R 2 + R 1 ]d x Rxdpt = 0, 



(48a) 
(48b) 
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FIG. 12: Model C with finite b = 0.5 for the same parameters as in figure [TT1 The plots are for times (a) t = 0, 3 and 6, and 
(b) t = 0, 3 and 7. For comparison, the corresponding profiles for b = are shown as dashed curves. Notice that the peak 
appears rather sharp due to the relative elongation of the vertical axis. Inset: For qualitative comparison, experimental results 
as given by Fig. 23 of Ref. ^3 for the surface velocity v and the local height h of the mobile layer which corresponds to R in 
our notation. A linear relation between v and R is observed, and the shape of R resembles that of our analytical result. (The 
direction of flow is inverted in the experiment compared to our model.) 



where the characteristic direction are the same as in the unperturbed case, = 1, (_ = — R\. Again all functions in 
the above equations have to be regarded as depending on a, /3. Derivatives with respect to x can be rewritten by the 
use of the relations 



d x /3 = - 



Ri 



l + Ri 



d x a 



ra(l + i?i) 



Using the latter relations, equation (|48a|) can be integrated with respect to a. The result is 



Z 2 (a,(3) 



da'- 



R' (a') 



_p l + R^a',0) 



(49) 



(50) 



By inserting this solution into equation (|48bfl one obtains an ordinary differential equation (with respect to (3) for 
i?2(a,/3). In formal analogy to equation (|34() its solution can be written as 



R 2 (a,P) = exp 



d/3' gi (a,/3') 



d/3'g2(a, /3') exp 



0' 



df3"gi(a,/3") 



where the functions g\ and g 2 are now given by 

rn 



92 («, P) 



(1 + Ri) [R Q {a) + m/R' {a) - Ri] 
1 



l + Ri 



R ig (a, 13) - d p ln(l + Ro(-/3)) + m 



da' 



R' (a>) 



(l + fli(a',/?))3 



(51) 

(52a) 
(52b) 



Most of the integrals in equation l|51() can be obtained in closed form. After a tedious calculation the final result can 
be expressed in terms of single integrals, that we give for completeness: 



R2(a,0) = 



1 + (-Ro(a) - Ri)R' (a)/m [ m 



R' (a) 



Ri — Ro(a) — In 



- a 1 (1 + i?! (a, /?'))(! + .Ro(-/?')) m 

R' Q {a')da' 



(l + R 1 )(l + R (-f3)) 
(l + i? (a)) 2 
R' (a') 



- [l + (R (a) + l)R' (a)/m 



/ a 'i + #iK/?) 

l + i?i(a,/3) l + Ri(a,-a') 



p (R 1 (a,(3) + l) 2 -(R 1 (a',P) + iy 



l + Ri(a',P) l + Ro(a') 



>(53) 
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The equations l|50|l and l|53(l provide the final result for the changes in the profiles due to a finite horizontal stress. 
The solution is valid for general initial profiles Rq(x). In figure IT21 we have plotted this solution for a Gaussian 
perturbation with the same parameters as in the previous section, cf figure 1111 

2. Physical discussion 

The main difference coming from finite b is the generation of a peak at the downhill front of the avalanche, both for 
positive and negative slopes m. From the model of Eqs. i|38|) one observes that the term proportional to b is controlled 
by the slope of the moving layer. Since this slope becomes steeper at the downhill front with evolving time, there is 
an increasing thickness of the mobile layer close to the front. From a physical point of view this amplification can be 
understood from the effect of the horizontal stress which increases with R. Thus at the downhill front there is a net 
force which pushes the material towards the front and induces an extra growth of the rolling layer. As demonstrated 
by the inset of Fig. 1111 our result is in qualitative agreement with measurements of the thickness of the mobile layer 
of grains along the symmetry axis of a triangular shaped avalanche moving on a static layer of limited thickness |ll| . 
For positive m there is, in addition to the homogeneous and linear decrease of Z ', a net transport of grains of the sand 
bed from the downhill front to the uphill end. 

V. CONCLUSION AND DISCUSSION 

In this paper, we have studied two sets of St Venant equations for the modeling of granular avalanches on an 
erodible bed. The models differ in the choice of the velocity profile within the flowing layer - either a plug flow (P) 
or a profile with a constant velocity gradient (£) — which give rise to different non-linearities. These models can 
been solved analytically by the method of characteristics, at least for sufficiently large stress anisotropy, and we have 
focused our attention on the situation where a uniform static slope is initially disturbed by a localized amount of 
rolling grains. We were able to compute the space and time evolution of the avalanche which either dies or grows, 
depending on whether the initial slope is below or above the angle of repose. 

Such unsteady and non-uniform flows are very demanding for the models as the description of fronts and the 
generation of shocks must be addressed. For the case of a plug flow we found that shocks occur at the uphill end of 
the avalanche above a positive critical value m c of the slope. Below this value, the asymptotic large time behavior can 
be computed, and we found that the amplitude and width of the avalanche grow linearly in time if m > 0, whereas 
for negative m the initial perturbation decays exponentially and the width of the deposited bump of grains scales like 
1/m. For a linear velocity profile, shocks occur for all slopes at a finite time that can be computed explicitly and is 
related to the shape of the initial distribution of rolling grains. By contrast to the previous case, shocks are located 
at the front of the avalanche. 

Our analysis shows that these models predict several interesting qualitative features of granular avalanches, that 
can be compared with experiments. However, in the present form they also certainly have a number of shortcoming. 
The plug flow assumption, for example, is not tenable for a rolling layer thickness which starts to be of the order of 
several grain diameters, as particles on the top of such a layer would not feel the damping due to the friction on the 
static bed. On the other hand, the linear velocity hypothesis yields a vanishing velocity U — > when the thickness 
of the mobile layer R — ► 0, which forbids any front to move. However, shock-less and well defined propagative fronts 
are observed experimentally in steady Q or unsteady |30| situations. Our model also does not correctly account for 
slope hysteresis which should differentiate between a starting and a stopping angle, s tart and ^ s top, respectively. In 
fact, the region between </> s top and </> s tart is precisely that of major interest as static slopes with an angle <fi > <p s top 
are stable but can generate growing avalanches when disturbed by a local amount of rolling grains. Because the 
non-linear terms in model V are all proportional to R, spontaneous avalanches (for R = 0) can never occur whatever 
the value of the initial slope. This would correspond to the maximal value of n/2 for </> s tart- This, however, implies 
that model V can be applied to the experimentally important regime of slopes which are only slightly larger than 
</>stop and sufficiently small compared to s tart- On the other hand, the uniformly growing solution for R in model C 
can be interpreted as evidence that both angles </> s tart and </> s top are identical so that beyond that angle avalanches 
occur even at R = (unstable slope). 

To include the above mentioned effects in the presented models, new directions have to be proposed. First, one 
can modify the dynamical equations itself while keeping their general structure (hyperbolic first order differential 
equations) with its analytic properties being tractable by the method of characteristics. Second, additional physical 
input could be used to either study the propagation of shocks beyond the shock time scale or one could implement 
boundary conditions on, e.g., the slope of the profiles to describe the dynamics close to the avalanche front. To be 
more specific, we discuss two points of particular interest: The velocity profile and hysteresis effects. 
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As proposed theoretically and demonstrated by experiments, an independence of the flow velocity must be kept to 
allow for thick avalanches. The commonly used strictly linear velocity profile with the rheological ansatz u(x, y) = jy 
(constant shear rate) leads to problems at the avalanche front since the depth averaged velocity U should stay finite 
and of the order of ^fgd when R — > 0. A possibility is to let 7 diverge in this limit. This would correspond to a 
crossover from model C to model V below a certain small value for R of the order of a grain diameter d (see |26J). 
The corresponding matching of characteristic curves is technically involved and moreover such a treatment would not 
provide an understanding of the physical mechanism of velocity selection. In fact, close photos of the foot of the 
avalanche fronts reported in 6] show a small gas-like region where grains are ejected from the dense flow. This means 
that this zone is, in a sense, outside of the present modeling framework and thus could allow for a discontinuity in 
R, for example, or for imposing an extra constraint on the profiles R and Z or their derivatives - e.g. a fixed slope 
of the free surface as observed on propagative fronts. Another more fundamental approach would be to consider the 
velocity U as an independent dynamical field which a priori is not related to R by a fixed velocity profile. Such kind 
of description has been applied to granular flow on a fixed plane 7] so that Z is not a dynamical quantity. 

Another feature of sand piles that should be included in the St Venant models discussed here is the hysteresis of 
avalanche dynamics as reflected by the existence of two critical angles. Experimentally it is observed that, at least for 
a fixed profile Z, the critical angles are actually functions of the thickness R of the mobile layer pfl. Since the friction 
coefficient \i is given by the tangent of the actual angle of the pile, a non-constant friction coefficient is expected. 
Indeed, at the scale of a single grain, the starting angle depends on the depth of "traps" due to the roughness of the 
static bed yjj ■ Therefore one should expect an increased value of \i below a typical thickness i?t rap (of the order of 
a grain diameter). This increased value determines then the staring angle </> s tart while the /i at large R corresponds 
to the stopping angle s t O p- A sufficiently large value of the friction coefficient at small R would lead, for the points 
where R < -Rt rap , to a "freezing" of the sand bed profile Z at its current value. 

Both modifications, constraints on the profiles at the avalanche front and a friction function fi(R), do not change 
the general structure of the dynamical equations studied in the present paper. Thus, the method of characteristics 
and our predictions should prove useful for a better understanding and modeling of non-stationary granular flow. 
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In this section a brief account of the general theory for systems of partial differential of first order and hyperbolic 
type is presented. Due to the relevance to granular flow problems we will concentrate on non-linear systems consisting 
of two equations for two functions of two independent variables. In the present context of this paper, the functions are 
R and Z and the i ndep endent variables correspond to space x and time t. For such systems a complete mathematical 
theory is available [3J. Our presentation will follow closely the latter reference. 

For hyperbolic systems the notion of characteristic curves is the central concept. Before introducing the general 
theory, we would like to motivate the introduction of characteristic curves or coordinates. This concept is particularly 
adapted to the case where the number of equations equals the number of independent variables. In the present case of 
two equations, the objective of the method of characteristics is to introduce instead of (t, x) a new coordinate frame 
(a, 0) so that along the two families of curves of constant coordinates a and (3 the partial differential equations reduce 
to ordinary differential equations with respect to a and f3. 

Let us demonstrate the method explicitly for the simple case of one linear partial differential equation for a function 
f(t, x) of the form 



The initial value will be prescribed at zero time, f(t = 0,x) = fo( x )- First, we have to define the characteristic curves 
[i(a),x(a)] where a is a variable parameterizing a given curve. These curves are specified in terms of their local 
tangent vectors (velocities), 
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APPENDIX A: METHOD OF CHARACTERISTIC CURVES 



a{t, x) d x f + b(t, x) dtf + c{t, x)f = 0. 



(Al) 




dx dt 



b{t, x) . 



(A2) 
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Integrating this equations yields a whole family of characteristic curves which is parameterized by the starting position 
x of the curves at t = so that x(a = 0) ~ Xq and t(a = 0) = 0. This ensures that from each position of the line 
t = exactly one characteristic curve originates, providing the trajectory along which the initial data fo(x) can be 
propagated in time. In order to see why the definition of Eq. (|A2|) is useful we compute the change of / along the 
characteristic curves, yielding 

£+c(i,x)/ = 0. (A3) 

The crucial observation is that the latter equation is just an ordinary differential equation, which is valid along the 
characteristic curves. To find the final solution f(t, x) to Eq. (|A1|) one proceeds as follows. First, one solves Eq. (|A2|) 
to obtain the relation between (t, x) and (a,a;o). Second, the ordinary differential Eq. I|A3(1 is solved with the initial 
condition f(a = 0) = fo(xo) which provides the solution f(a,x ). Finally, the parameters a and Xq are computed for 
a given coordinate (i, x) to get the solution /(£, x) in the original coordinate frame. 

Having outlined the general idea behind the method of characteristics, we can go ahead and turn to the case of two 
non-linear hyperbolic equations for two functions. We will consider a general system of the form 

Li = A 1 d x R + B 1 d t R + C 1 d x Z + D 1 d t Z + E 1 =0 (A4a) 
L 2 = A 2 d x R + B 2 d t R + C 2 d x Z + D 2 d t Z + E 2 = (A4b) 

for the functions R(t, x) and Z(t, x), where the coefficients A±, A 2 , B±, . . ., are known functions of x, t, R and Z. The 
type of this system depends on the coefficients. For the hyperbolic case in which we are interested here one needs 
that ac — b 2 < with the functions 



a=[A,C], 2b = [A, D] + [B, C], c=[B,D], (A5) 

where [X, Y] =X 1 Y 2 - X 2 Y X . 

The goal is again to reduce the above system to a system of ordinary differential equations with respect to new 
coordinates a, [3. Since we have now to deal with two unknown functions R(t, x) and Z(t, x) we start by searching for 
a linear combination L = X\Li + X 2 L 2 of the differential operators in Eq. (IA4|) so that the derivatives of R and those 
of Z combine to derivatives in the same direction. These directions will be the velocity vectors of the characteristic 
curves and thus determine the new coordinate frame. Let us represent an arbitrary curve in the x — t plane by 
[t(o~), x(a)] with a denoting the running parameter along the curve — note that a finally will play the role of a or (3. 
Then the condition that in L both functions R and Z are differentiated in the tangential direction of this curve reads 

dx/da _ XiAi + X 2 A 2 _ AjA + A 2 C 2 (M) 
dt/da ~ Aifli + X 2 B 2 ~ AiDi + X 2 D 2 ' 1 ' 

Next, we consider the change of the functions R and Z along the curve [t(a), x{a)\. It is given by dR/da = d x Rdx/da+ 

dtRdt/da and analogous for Z. Multiplying L with either dx/da or dt/da and using the conditions of Eq. I|A6(I one 
gets 

= (X 1 A 1 + X 2 A 2 ) ^ + (Aid + A 2 C 2 ) ^ + {\ l E l + X 2 E 2 ) ^ (A7a) 
da da da da 

-j-L = (X 1 B 1 +X 2 B 2 )^ + (X 1 D 1 + X 2 D 2 )^ + (X 1 E 1 +X 2 E 2 )^. (A7b) 
da da da da 



If the functions R and Z satisfy the system of differential Eqs. (|A4|) we have L = 0, and we obtain the following four 
homogeneous linear equations for the coefficients Ai and A 2 which result from Eq. I)A6|I and Eq. 1A7(I . 
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This system is obviously over-determined. Thus, in order to have a non-trivial solution, the determinant of every pair 
of rows in the matrix of coefficients of Ai and A2 has to vanish. The relations following from this conditions are called 
characteristic relations. 
In particular, from the first two Eqs. IjASajl . (| A8 b|) . one obtains the condition 



dt\ 2 „ , dx alt ( dx 



da ) da da \ da ) ^ ^ 

with the coefficients given by Eq. (| A5|) . From this condition it becomes clear why the method of characteristics 
applies to hyperbolic systems. For those systems we have, as mentioned above, ac — b 2 < 0, and Eg. (|A9fl has two 
different solutions and thus two different characteristic directions (dx/da, dt/da) through each point. In the following 
we assume, without any restrictions, that a ^ so that dx/da 7^ and we can introduce the slope 

The two different real solutions of aC/ 2 — 2b( + c = for these so-called characteristic directions will be denoted by 
and £_ , respectively. These characteristic directions are in general functions of t, x, R and Z. Two one-parameter 
families of characteristic curves follow from the directions by integration of the ordinary differential equations dt/dx = 
£+(x,t, R, Z) and dt/dx — £_ (x, t, R, Z). In the following we will denote the families of curves by C+ and C_. These 
two families of curves define a curved coordinate net if the curves arc represented as a(x, t) = constant and (3{x, t) = 
constant for family C_ and C+, respectively. The functions a(x,t) and (3{x,t) are called characteristic parameters. 
The coordinates (t, x) corresponding to a given pair (a, 0) can be obtained as follows. Consider a curve X given by 
[x(s),i(s)] which has nowhere a characteristic direction as tangential vector. In practice, X will be usually the line 
t = where the initial data are defined. In addition boundary conditions, e.g., a fixed flux at a given position, can 
be specified by a curve with x ^constant. Through the two points s — a and s = (3 on the curve X one follows 
the characteristic curve of family C_ and C+, respectively, up to the point where the two curves intersect. The 
new coordinates of this intersection point (t,x) are then (a,f3). The characteristic parameters a and (3 can now 
replace the parameter a for the curves of family C+ and C_, respectively, so that one has dt/da = C + dx/da and 
dt/d/3 = C-dx/dfi. 

Next, we have to find equations which determine the evolution of the functions R and Z along the characteristic 
curves. This can be done by eliminating Ai and A2 from the Eas. l|A8a|) and i|A8c|) . Using dt/da = (dx/da, where £ 
denotes either C+ or £_ and a is either a or /?, one obtains 

T-j- + K - S)^ + (KC -H)^=0, (All) 



with the coefficients 



da da da 



T=[A,B], S=[B,C], K = [A,E), H=[B,E]. (A12) 



If we apply the latter equation to the curves of C+ and C_ and combine them with the equations for the characteristic 
curves, we finally obtain the following four characteristic equations which are differential equations for the four 
functions x(a,/3), t(a,/3), R(a,(3) and Z(a,(3) and replace the original system of Eq. (|A4(I . 

d a t — (+d a x = , (Al3a) 

dpt-t-dpx = 0, (A13b) 

Td a R + {a(+ - S)d a Z + {K( + - H)d a x = , (A13c) 

TdpR + (aC- - S)d f3 Z + {KC- - H)d p x = . (A13d) 

All the coefficients in this system are known functions of x, t, R and Z. It can be shown that every solution of this 
characteristic systems satisfies the original system of Eq. (|A4(1 provided that d a xdpt — dpxd a t — (£_ — C, + )d a xdpx is 
non-zero. With the derivation of Eq. I|A13|) we reached our initial objective to reduce the partial differential equations 
to a form which resembles that of ordinary differential equations along certain curves. This can be seen from the fact 
that each equation contains derivatives with respect to only one of the coordinates a and f3. Moreover, the system 
has to the convenient property that the coefficients do not dependent on the independent variables a and f3. 

Now we are in the position to outline the strategy for solving an initial value problem for the system of Eq. i|A4fl . 
Let us assume that the initial values of the functions R and Z are given on the line t = by Ro(x) and Zq(x), and 
that this line has no characteristic directions. This line corresponds then to the curve X introduced above. We may 
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consider this curve as the image of the characteristic parameters obeying the relation a + (3 = 0. Then we have to 
solve the system of Eq. (| A13I) with the initial conditions 

t(a, — a) — 0, x(a, — a) = a, R(a,—a) = Ro(a), Z(a,—a)=Zo(a). (A14) 

Due to the particular simple structure of the system of Eq. (|A13|) . this problem can be treated as completely as 
the initial value problem for ordinary differential equations. It is this formulation of the non-linear hyperbolic flow 
problem which we used throughout the paper to solve the partial differential equations exactly. Finally, we note that 
this method can be generalized to an arbitrary number n of equations. Then the system will have n characteristic 
directions £„ and correspondingly n different families of characteristic curves. However, the n resulting characteristic 
parameters can no longer be interpreted as a new coordinate frame since there are only the two coordinates t and x. 



APPENDIX B: DERIVATION OF THE SHOCK CONDITION 



In this appendix we first review the mechanism for the generation of shocks and their mathematical definition. 
Then we provide for model V the details of the calculations for the shock existence criterion and the time and position 
of the shock. For model C shocks are always generated, and we derive the simple result of Eq. I|45|l for the shock time. 
In Appendix [3] we have assumed that the characteristic curves of one family (either C + or C_) do not intersect. 
Only if this is true there is a well defined mapping between the original coordinates (t, x) and the characteristic 
parameters (a, (3). However, depending on the initial data at zero time, it is possible that characteristics of the same 
family intersect at a finite time. Beyond this shock time the system of partial differential equations fails to have a 
single valued solution but only multi-valued solutions or even no solution at all exists at later times. The points of 
intersection of characteristic curves are enclosed by an envelope, cf. Fig[SJ The earliest time where a shock appears is 
the position of the cusp of this envelope. Technically, the envelope is defined by the condition that for every position 
on the envelope there exists a characteristic curve that touches the envelope at the position so that both curves have 
the same tangential direction. If we represent the envelope as [t e (a), x e (a)] where a is used as the parameter changing 
along the envelope then the conditions read 



x a (t e (a)) = x e (a), d a x a {t)\ t=ta{a) = 0. 



(Bl) 



where x a (t) is the trajectory of the characteristic curve along which a is constant. The second condition follows from 
the requirement that the tangent vector [1, dtx a (t)] of the curve x a (t) is parallel to the tangent [dt e (a) / da, dx e (a) /da] 
of the envelope. To see this, one takes the derivative of the first condition of Eq. IjBlfl with respect to a so that one 
obtains 



dx e (a) ,,,dt e (a) 

= Vtx a {t)—j— + d a x a (t)\ t=te{a) 



da 



da 



(B2) 



The collinearity of the two tangent vectors requires then the last term on the rhs of Eq. (|B2|I to vanish. Notice that 
in model V the characteristic curves along which [3 is constant are straight lines, and thus they can never produce a 
shock. 

In the following we focus on the initial profile Rq(x) given by Eq. @ which allows for an explicit calculation of the 
envelope and the condition for shocks. For this profile it is useful to consider three different sectors in the t — x plane, 
see Fig. 0] Using the above conditions, one obtains for the envelope in sector (II) the result 



t e (a) 
x e (a) 



1 



m — 1/5 
1 

m — 1/5 



In 



1 



h'(a) h(a) 
h(a)e {m ~ 1/S)t ' {a) +ln 







+ ma — ro 







(B3a) 
(B3b) 



where the function h{a) is given by Eq. (|18fl . In sectors (I) and (III) the conditions of Eq. (|B1|) cannot be satisfied 
and thus the characteristics in these sectors never form an envelop. This can be seen as follows. In sector (I), the 
characteristics along which a is constant are given by 



x a {t) = a + Ro(a) 



e (m+l/S)t _ ^ 



m + 1/5 

Using this expression, the second condition of Eq. I|B1(I becomes 

R' (a) 



= 1 



m + 1/5 



,(m+l/5)t 



(B4) 



(B5) 
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For to + 1/5 > 0, the term added to one on the rhs is positive for t > since Rq(o) > for negative a in sector (I). In 
the opposite case to + 1/(5 < 0, the same argument applies since the expression in the square brackets in now negative 
for t > 0. This shows that the rhs is always larger than one, and the condition is never fulfilled. In sector (III) we use 
a different argument to show that no a=constant characteristics, which originate from positive x at t — 0, form an 
envelope. From the characteristics in sector (III) , the conditions of Eq. (|B1|I are formally fulfilled by the expression 



t e (a) 
x e (a) 



1 



In 



TO — 1/(5 

_ Ro(a) 
a KM 



1 



1/6 



R' (a) 



(B6a) 
(B6b) 



for the envelope. However, it remains to be checked that this curve is indeed located in sector (III), i.e., if its 
coordinates are larger than the boundary between sectors (II) and (III), cf. Eq. (|14b(l . which yields 



x e (a) > 



m 



1/(5 



,(m-l/S)t e (a) 



(B7) 



Using the definition of Rq(x), cf. Eq. 
turns out to be equivalent to 



©, and the relation R' (a) = — (l/6)Ro(a)/(l + Ro(a)), the latter condition 
■R {aY 



In 



Ro{a) 



< 1. 



(B8) 



Since Ro(a)/ro < 1 for a 7^ 0, this condition is in fact never fulfilled which proves the absence of shocks in sector 

nil). 

Knowing that shocks, i.e., the cusp of the envelope, can occur only in sector (II) we can try to obtain the condition 
for shock generation and the time and position of the shock. First consider a negative slope to < 0. Then the 
characteristics in sector (II) saturate at large times, 



lim x a (t) 



ro 



1/6- m 



1 



-a 



In 



ro 



ro 



h(a) 



(B9) 



Since h(a) is a monotonously decreasing function for negative a, the expression in the square brackets is monotonously 
increasing in a. Thus the characteristics retain the original order for all times, i.e., they never intersect. The situation 
is more complicated for positive to. Let us assume that there exists a finite value m c so that only for to > m c shocks 
are produced. Then, one expects that at to = m c the time t s for the shock appearance tends to infinity in order to 
have no shocks at finite times for to < m c . Since the shock position (t s ,x s ) is the cusp of the envelope, we have to 
analyze the large time behavior of the envelope of Eq. (|B3(I . We start with the assumption that the critical value 
to c < 1/(5, and, in fact, at the end we will confirm this assumption. For to < 1/6 one has the asymptotic behavior 
t e {pt) ~ —a, and thus we consider large negative values for a in the following. The shock time is given by the minimal 
time of the envelope t s = t e (a m ) with a m the parameter at the minimum, i.e., dt e (a) / 'da = for a = a m . Then 
close to the critical slope to c , we expect that a m — + —00. For large negative a, the function h(a) of Eq. (|18fl has the 
asymptotic form 

h(a) ~ rl^i-maye^ 1 -^^/^ 

(|B3|I the condition dt e (a)/da 



with v — (2/(5)/ (to + 1/(5). Using this expansion in Eq. 
large a independent of a and assumes the simple form 



(B10) 

becomes at asymptotically 



v- 1 



(Bll) 



Since we assumed that the shock time t s — > 00, this condition has to be regarded as an equation for the critical slope 
m = m c . Since v depends on to the equation is quadratic in to, and it has one negative solution which we have to 
discard and the other solution gives the critical slope beyond which shocks occur, 



V2-1 



(B12) 



This is the result given in Eq. Ijl9|l . The behavior of t s close to to = m c can be obtained by computing the leading 
correction to the (constant) asymptotic expression for dt e {a) / da. We find a correction *~ 1/a which in turn yields 
the leading order of a m close to to c , 
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Since at large negative a the time coordinate of the envelope behaves as t e (a) ~ —a and t s = t e (a m ), we obtain the 
following power law for the shock time close to criticality, 

t s , (B14) 

to — m c 

as given by Eq. H2U|) . The precise time t s and position x s of the shock is given for to sufficiently close to m c by the 
envelope of Eq. (|B3() at a = a m of Eq. i|B13(l . At larger to the coordinates (t s ,x s ) are difficult to compute. However, 
at sufficiently large to closed formulas can be obtained. The reason for this is that for to larger than some threshold 
the minimum of t e {a) is always at a m = 0, i.e., the shock is located on the boundary between sectors (II) and (III). 
To see this, we expand the envelope of Eq. I|B3|) now around small negative a. This gives 



f \ 1 fmS(l+ro) — 1\ 1 — 2r <5m 2 

*e(« = 777 M " + 77T VA~TTI S 1 x a + ° a B15a 

to -1/(5 V r / (1 + r )(om(l + r ) - 1) 

1 \ 1 — 2r <5m ^ 2 



xJa) = r 8[l + —)+ — a + 0(a 2 ) . (B15b) 

\ r J 1 + r 

The minimum of t e (a) is at a = if the coefficient of a in Eq. (|B15a|l is negative. The denominator of this coefficient 
has to be positive since otherwise the argument of the logarithm in Eq. i|B15a|l would be negative. Thus if the slope 
to fulfills the two conditions to > l/(2ro#) and m > 1/(5(1 + ro)) simultaneously then the shock position is given by 
(t s ,x s ) = (t e (a = 0), x e (a = 0)). As mentioned already in Sec. IIII Al for ro < 1 the first condition is relevant whereas 
for ro > 1 the latter condition dominates. Now one may ask if it is possible that m c is larger than latter thresholds 
so that the shock would occur no longer at a m of Eq. l|B13p but at a m — 0, i.e., on the boundary between sectors (II) 
and (III) at rather small times. In fact, for ro < 1 the condition m c > 1/(2tq5) is never fulfilled so that the minimum 
remains at a m of Eq. I|B13(I . For ro > 1 the condition m c > 1/(5(1 + ro)) leads to ro > v2. In the latter case the 
shock occurs always at a m = and the new critical slope is given by 1/(5(1 + ro)). However, it should kept in mind 
that the width of the perturbation Rq(x) of Eq. @ is proportional to 5 only for ro < 1. For larger ro > 1 the width 
is proportional to ro<5. From the first term x s =x e (a = 0) in Eq. (jB15bjl we thus conclude that the shock occurs at 
the uphill end of the avalanche with the shock position approximately given by the uphill end of the perturbation at 
t = 0. 

Finally, we study the generation of shocks for model L. We do this by using an approach which is more adapted 
to the special structure of the solution of this model. We do not use directly the definition of Eq. QBlfl but look for a 
discontinuity in the profile R(t, x) as a function of x. If there is a jump in R(t, x) at some position x then a shock is 
generated and the earliest time where this happens if the shock time t s . We start from Eq. (|42|l which gives 

R(a,0) = Ro(a) + mt. (B16) 

By taking the derivative with respect to x, one obtains 

d x R=R' (a)d x a. (B17) 

Since R'q(x) remains finite, we have to search for a divergence in d x a. The characteristic parameter a(t,x) can be 
obtained from Eqs. I|41c|) and H42|) which yield 



TYi 

a(t,x) = x + —t 2 + Ra(a)t, (B18) 



which leads to 



1 — R (a)t 

Since there is always an interval of values for a where a localized Ro(x) has a positive derivative, shocks are generated 
always. The time scale t s for the occurrence of the shock is the earliest time where d x a diverges, i.e., it is given by 
the maximal slope of the initial perturbation, 

*- = — ' ( B2 °) 

maxK (a;) 

which is Eq. I|45|) . The position x s of the shock follows from Eq. (|B18Jl . For the Gaussian perturbation Ro(x) = 
ro exp(— x 2 /5 2 ) discussed in Section HV Al one has t s — y^e/25/ro and the position is given by 



e 

4 \r 



2 



= -V26- - A [ — ) to. (B21) 
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Note that the shock occurs at the downhill front of the avalanche as opposed to the uphill position in model V . 
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